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Abstract 

A Markov process is constructed to numerically study the phase sepa- 
ration in the 6-vertex model with domain wall boundary conditions. It is a 
random walk on the graph where vertices are states and edges are elemen- 
tary moves. It converges to the Gibbs measure of the 6-vertex model. Our 
results show clearly that a droplet of c vertices is created when Boltzamnn 
weights are in the antisegnetoelectric region. The droplet is a diamond- 
like shaped curve with four cusps. 

Introduction 

It is well known that the 6-vertex model is exactly solvable and has phase 
transitions. The history and the classification of phases in the 6-vertex model 
as well as many interesting facts about the structure of the partition function 
of the 6-vertex model with periodic boundary conditions can be found in |10|. 

m 

There is an important function of Boltzmann weights of the model which is 
usually denoted by A (see [SI; CO] an d section[3J). The 6-vertex model with peri- 
odic boundary conditions has 3 phases in the thermodynamical limit, depending 
on the value of A. One is the totally ordered (frozen) phase, with A > 1, the 
second is the disordered (critical) phase, with |A| < 1, and the third is the 
partially ordered (antisegnetoelectic) phase with A < — 1. 

The 6-vertex model with domain wall boundary conditions on a square Nx N 
grid perhaps was first considered in [?] in the process of computation of norms 
of Bethe vectors. The partiction function of this system can be written as the 
determinant of a certain N x N matrix Its asymptotics in the thermody- 
namic limit N —> oo were analyzed in \7\. It is related to matrix models, which 
was pointed out and exploited in |17| . 

The 6-vertex model with A = 1/2 is also known as the ice-model. This model 
with domain wall (DW) boundary condition is closely related to the enumeration 
of alternate sign matrices |§J . It also has other interesting combinatorial features 
(see for example (TB) . [T5]). When A = the 6-vertex model is equivalent to the 
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problem of counting of weighted tilings of the Aztec diamond (see for example 
|7], |Hj and references therein). 

The spatial coexistence theory of different phases and the interfaces sepa- 
rating phases is an important part of statistical mechanics. Growth of crystals 
is one of the well known phenomena of this type. This is also closely related to 
the limit shape effect in statistics of Young diagrams an d plane partitions. 

In dimer models related to enumeration of plane partitions and domino 
tilings, the interface between the disordered and totally ordered phases is also 
known as an arctic circle phenomenon p^. 

Dimer models on bipartite planar graphs with periodic weights are exactly 
solvable models where this phenomenon has been studied in '11 . In dimer 
models the limit shapes or interfaces (curves separating phases), under broad 
conditions, are real algebraic curves Since at A = the 6- vertex model is 
equivalent to a dimer model these results imply that such a phenomenon exists 
in the 6-vertex model for A = 0. The natural question is whether the spatial 
coexistence of phases happens only at the free fermionic point or if it occurs for 
all values of A. Numerical evidence suggesting the existence of a limit shape in 
the 6-vertex model with domain wall boundary conditions for all weights was 
obtained in [Tl] . 

Here we report results of numerical study of the 6-vertex model with DW 
boundary conditions in all phases of the model. Our method is different from 
|14|. To generate a random configuration in the 6-vertex model we construct the 
Markov process which is equivalent to a random weighted walk on the graph 
where the vertices are states of the model and edges are local moves which 
transform states into other states. This process satisfies the detailed balance 
condition and therefore converges to the Gibbs state of the 6-vertex model. It 
is also known as Monte-Carlo with local update and as a heat-bath algorithm. 
In statistical mechanics such processes are known as Kawasaki, or Glauber dy- 
namics. For a more effective version of this algorithm known as the "coupling 
from the past" algorithm, see |12| . 

For periodical boundary conditions the system can be either in the ordered 
(segnetoelectric) phase, disordered (critical) phase, or antisegnetoelectric (non- 
critical) phase, depending on the values of Boltzmann weights. 

Our results confirm the conclusion from ^ and \TQ that there is a coexis- 
tence of ordered and disordered phases in the 6-vertex model. They also clearly 
indicate that for A < —1 there is a coexistence of all three phases. The outer 
layer is an ordered phase. It follows by the ring of disordered phase. Finally, 
there is an inner droplet of the antisegnetoectric phase. This phenomenon was 
first conjectured in ^1] using a different numerical method. The shape of the 
inner droplet has four cusps and is reminiscent of one of the limit shapes for 
dimcrs on a square-octagon grid equivalent to the diablo tiling . 

Acknowledgments. We are grateful to R. Kenyon and A. Okounkov for 
many illuminating discussions, to K. Palamarchuk for valuable comments, and 
to T. Yates for help with the implementation of the algorithm in the C language. 
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1 Weights and local moves 



1.1 States 

States of the 6-vertex model on a square lattice are configurations of arrows 
assigned to each edge. The 6-vertex rule is that the total number of arrows 
coming into any vertex should be equal to the total number of arrows going out 
of this vertex. Each configuration of arrows can be equivalently regarded as a 
configuration of empty edges (arrows oriented South-North and East- West) and 
occupied edges, or thick edges (arrows pointing in the opposite directions). It is 
clear that thick edges will form paths. Possible configurations of paths around 
a vertex are shown on fig. 11.11 
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Figure 1: The 6 vertices and their weights 



We will use oi, <Z2, b\, 62, C2, and C2 as names of the vertices. We denote by 
the same letters Boltzmann weights assign to these vertices. 
Domain wall boundary conditions are shown on fig. [21and|3J 
For domain wall boundary conditions every path in a 6-vertex configuration 
will have one end at the North boundary of the square and the other end at the 
East boundary of the square. These paths can be regarded as level curves of a 
height function. The lowest height function is shown on fig. Eland the highest 
height function is shown on fig. [21 

1.2 Weights 

The weight of a state is the product of weights of vertices and the weight of a 
vertex is determined by rules from fig. 11.11 The partition function is sum of 
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Figure 3: Domain Wall Boundary Condition Low 
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weights of all configurations: 



Zq v = y. I I w (vertex, state) 



states vertices 



where w(vertex, state) is the weight of the vertex (see fig ll.f \ . 
The ratio 




is the probability of the state. This is the Gibbs measure of the 6- vertex model. 
1.3 Local moves and the graph of states 

Now let us describe local moves in the space of states. Such a move changes the 
configuration of arrows at the minimal number of edges near a given vertex and 
it acts transitively, i.e. any given state of the model can be transformed to any 
other given state of the model by a sequence of such moves. 

Such moves are most transparent in terms of height functions. There are 
two types of local moves: 

• The path from fig. Q]we can move up, i.e. to the path from fig. El We 
will call this move flip up. 

• The path from fig. 0can be moved down, i.e. to the path from fig. 0] We 
will say that this is the flip down. 



Figure 5: The local configuration after an up flip has occurred 

Such moves with all possible surrounding configurations we will call flips up 
and down. 

For each flippable vertex we introduce effective weight as follows: 




Figure 4: A local configuration that may be flipped up 
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• For a vertex flippable up the effective weight is the product of weights of 
all vertices that can be affected by the flip, i. e. the vertex itself, and the 
neighboring vertices to the North , the East, and the North-East of it. 

• Similarly for a vertex flippable down the effective weight of it is the product 
of weights of the vertex itself, and of next neighboring vertices to the 
South, West, and South- West of it. 

The effective weight is always the product of four factors. The effective 
weight of vertex v in the configuration S we denote by W V (S). 

2 The Markov process 

2.1 General strategy 

Consider the abstract graph with vertices being states of the model and with 
edges being local moves. This graph is clearly connected. Our goal is to con- 
struct a random walk on this graph converging to the probabilistic measure 
vertices of this graph which is the Gibbs measure of the 6-vertex model with 
DW boundary conditions. 

Let us recall some basic facts. Let T be a finite connected graph and q : 
V(T) — > M.+ be a probabilistic measure on the set of vertices of T. Let M = 
{p{a — > &)} a ,fcev(r) be the matrix of the Markov process describing a random 
walk on T. A traveller moves from a to 6 with the probability p(a — > b). 

The matrix M must satisfy the total probability condition: 

5>(a^6) = l 

b 

If in addition it satisfies the detailed balance condition 

q(a)p(a -> b) = q(b)p(b -> a) 

then it is known that the Markov process converges to q. For details about 
Markov sampling and estimating convergence times, see |13j . 

Now our goal is to construct such random walk converging to the Gibbs state 
of the 6-vertex model. At some point the rate of convergence of this Markov 
process becomes an important issue. To avoid the complicated analysis of es- 
timating mixing times we will modify the algorithm and will use the "coupling 
from the past" version. This will be explained later. 

2.2 The Markov process for the 6-vertex model 

We want to construct Markov process which chooses a vertex at random, 
then with the probability which we will describe below it will either flip the 
configuration up at this vertex, or will flip it down, or will do nothing. The 
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probability of passing from the state S a to the state Sb in this process can be 
w written as follows: 

P(s a ^Sb) = w ^-Y p ^^ s ^ W 

V 

= ^-1 E s{s a ,s b)+ £ Ms a ^s b)+ 

\ (v)non-flip (v) flip— up 

(2) 

]T Pv(S a =^Sb)+ Y Pv(S a =>Sb)\ 

(v) flip— down (v)bi — flip / 

only 



\. — f (# non-flippable mS a )S Sa .s b + V P v (S a =► S b )+ 
# vertices \ ^— ' 

\ (v) flip— up 

only 

flip — down (v)bi— flip / 



(v) flip — down (v)bi— flip 

only 



#non- flippable ^ ^ 
# vertices Sa,Sb 

#flippable / 1 



Pv{S a => Sb)~ 



#verticcs \ #flippable 

\ (v) flip— up 

only 

—1— ]T p„(s a =»s„) + 1 £ W=^ s *) 

#flippable ^ #flippable ^ 

(v) flip— down {v)oi—flip / 

only 

(4) 

Here #flippable is the number of flippable vertices and #verices is the total 
number of vertices. 

Algorithmically, this means that we do the following: 

1. With probability P = #n °"; fl r ;PP a s blc , do nothing (that is, restart the loop.) 

2. With probability P = # fli P-"P-°"'y + *^^ s n ~ only + #bi - flip = gg^ , 
continue to the next part. 

If the algorithm continues, select a flippable vertex with the probability: 

P ^ Mion ) = #flip P ables (5) 
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At this selected vertex the configuration can be either flippable only up, or 
only down, or in both directions. Depending on this proceed according to the 
following rules: 

Three possible conditions now exist: 

1. The vertex is flippable down only. Two options: 

• Flip vertex down with probability P v (S a Sb) = pW v (Sb) 

• Stay with probability P v (stay) = 1 — pW v (Sb) 

2. The vertex is flippable up only. Two options: 

• Flip vertex up with probability P v (S a ==>■ Sb) = pW v (Sb) 

• Stay with probability P v (stay) = 1 — pW v (Sb) 

3. The vertex is flippable up and down. Three options: 

• Flip vertex down with probability P v (S a => Sb) = pW v (Sb) 

• Flip vertex up with probability P v (S a =>■ S b >) = pW v (S b >) 

• Stay with probability P v (stay) = 1 - pW v (Sb) — pW v (S b i) 

Here W v (S b >) W v (Sb) are the effective weights of the vertex v in the states 
obtained by flipping up or down at this vertex from the state S a . Effective 
weights were described in section fl. 31 

The parameter p is chosen such that all probabilities of transitions should 
be positive. In other words it should satisfy all conditions p < w ^ s ,^ where v is 
a vertex flippable in the state S either up or down, but not biflippable and S' 
is the configuration after the flip. At every biflippable vertex in the state S we 
should have p < w (s>)+w(s") wri ere S' is the result of the flipping S up at v 
and S" is the result of the flip down. 

This process satisfies the detailed balance condition, and the total probability 
condition. Since the graph of states with edges being local moves is connected, 
this process converges to the Gibbs state of the 6-vertex model. The process 
also depends on the choice of p. It slows down when p is small. 
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3 Random states in the 6-vertex model with 
DW boundary conditions 



3.1 Phases in the 6-vertex model 



3.1.1 



One can write weights of the 6-vertex model as 
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(6) 



(7) 



(8) 



Here Ei,E2, and £3 are energies of the interaction of arrows (or energies asso- 
ciated with the local shape of level curves of the height function) and H x and 
H y are magnetic fields. 

In this interpretation arrows can be regarded as spins interacting with the 
magnetic filed such that the energy of a vertical arrow is ±H X depending on 
whether the arrow is heading up or down. The energy of a horizontal arrow is 
±H y depending on whether it is oriented left or right. We assigned the energy 
of an arrow to the energies of adjacent vertices. 

Notice that since the total number of c\- and C2-vertices satisfy the relation 
n(c\) — n{c%) = N the partition function changes by an overall factor only when 
we change c\/c 2 . The total numbers of a and b vertices satisfy similar relations: 
n{a\) = 71(02) and n{b\) — n(b 2 ). Because of this for the square lattice with 
DW boundary conditions we can set a\ = a 2 = a, b\ = b 2 = b, and c\ = c 2 = c 
without loosing generality. 



Let us recall the phase diagram of the 6-vertex model ^UJi EH with periodic 
boundary conditions in the absence of magnetic fields. The important charac- 
teristic of the model is the parameter 



The phase diagram for the 6-vertex model with periodic boundary conditions 
in the absence of magnetic fields is shown on fig. |S| 
There are four phases: 

1. Phase I: a > b + c(A > 1). This an ordered phase where there are two 
possibilities for the ground state. It either consists of ai-vertices or of 02- 
vertices. In either case any change in the ground state gives the state with 
the total number of b and c vertices comparable with the linear size N of 
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Figure 6: The phase diagram for the six- vertex model in terms of the weights 
of a and b, assuming c = 1 . 



the system. Thus, as N — ► oo the energy of these two ground states is 
macroscopically separated from the energy of other states. In other words 
these are two frozen ground states. 

2. Phase II: b > a + c (A > 1). This is an ordered phase with double 
degeneracy of the ground state. The first possibility is when all vertices 
are b\ vertices, the second possibility is when all vertices are 62-vertices. 
As in case of phase I, this is a frozen phase. 

3. Phase III: a,b,c < 2±|±£(|A| < 1). This is a disordered phase. Local 
correlation functions decay as a power of the distance in this phase. These 
are the values of a, b, c when | A| < 1. In particular, the free fermionic curve 
A = lies entirely in this phase. It is shown by the dotted segment of the 
circle on fig. EJ 

4. Phase IV: c > a + b (A < — 1). This is an ordered phase with so-called 
antisegnetoellectic ordering (see fig. Q). The ground state in this case 
consists of alternating c\ and c 2 vertices. It is double degenerate due 
to the breaking of Z 2 -translational symmetry. In this case microscopic 
deviations from the ground state are possible. There is a finite correlation 
length in the system and local correlation functions decay exponentially. 

For details about phase transitions, magnetization, and the antiferroclcctric 
phase etc. see and [3]. 

3.2 The structure of a random state 
3.2.1 Free fermionic point 

This is the case when A = 0. It is convenient to parameterize weights in this 
case as 

a = pcosu, b = p sin it, c = p. 

When a = b = 1/ V2 this model is equivalent to the domino tiling of the Aztec 
diamond. The limit shape was computed analytically in and is a circle. 
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Figure 7: Antisegnetoelectric phase - in this phase, zig-zag paths alternate with 
zig-zags formed by empty edges 



The height functions of the average states for several values of the parameter 
u are shown on fig. [S] For A = 0, the limit shapes can be computed explicitly 
using methods of |11| and they are ellipses, which agrees with fig. [SJ 




2a = b a = b a = 2b 



Figure 8: Free fermionic point with A = 



3.2.2 Ordered phases 

In phase I the a-vertices dominate and the Gibbs state in this case is given by 
the lowest height function fig. |31 

In phase II the b vertices dominate and according to 1 1 Oj we should expect 
that the average state will be the state with the domination of 6-vertices. In 
other words the average state in this case is given by the highest height function 
ng.E 

3.2.3 Disordered phase 

In this case it is convenient to use the following parametrization of weights: 
a = rsin(7~ u), 6 = rsinu, c = rsin7 
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with 0<7<7r, 0<u<7, and r > 0. In this parametrization A — cos 7. 

Phase III contains the free fermionic curve A = 0. Since all this phase is 
critical one may expect that the nature of the Gibbs states will be similar for 
all parameters a, 6, c in this region. In particular, one can expect the existence 
of the limit shape as in the case A = 0. The particular form of the limit shape 
may vary but the following common features should common for all values of 
a, 6, c in this region: 

• The limit shape is a smooth curve having exactly one common point with 
each side of the square. At this point the limit shape is tangent to the 
side of the square. 

• Inside of the boundary of the limit shape the height function is a smooth 
function and it has continuous first derivative at the boundary. The second 
derivative has a discontinuity in the normal direction to the boundary of 
the limit shape. 

• Outside of the boundary of the limit shape the height function is linear. 
Examples of Gibbs states in the disordered phase are shown on fig. HO ^| 

EH 
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2a = b a = b 

Figure 9: Disordered phase with 7 = j 



a = 2b 



3.2.4 The antiferroelectric phase 

This region c > a + b is the one which is non-critical and which is also not 
ordered. In the periodic case the ground state has the domination of c-vertices 
as it is shown on fig. IT21 

It is convenient to use the parameterization 

a = r sinh(?7 — u) , b = r sinh u, c — r sinh r\ 
with < u < 77. In this parameterization A = — cosh 77. 



12 



2a = b a = b a = 2b 

Figure 10: Disordered phase with 7 = § 




2a = b a = b a = 2b 

Figure 11: Disordered phase with 7 = § 



In the case of DW boundary conditions there is a competition between very 
rigid restrictions on the states near the boundary which allows only a and b 
vertices near the boundary and the tendency of the system to have as much as 
possible of c vertices. 

Numerical simulations show that these competing tendencies resolve in the 
separation of three phases. It is fairly convincing from the fig. ^] that the 
following should take place: 

• The system forms a macroscopical droplet of the antiferroelectric phase 
with the boundary that does not touch the square. The height function in 
this domain is linear. The boundary of this domain has four cusps pointing 
towards sides of the square lattice. This phase is noncritical. Correlation 
functions in this region decay exponentially. 

• Near the boundary the system is ordered. This ordered region is bordered 
by the disordered region where the height function is smooth. The disor- 
dered phase is critical. There is a finite magnetization, which means there 
are excitations with linear dispersions and therefore correlation functions 
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decay according to a power law. The boundary between ordered and dis- 
ordered phases is a smooth curve with the features similar to the |A| < 1 
case. 




2a = b a = b a = 2b 



Figure 12: Antiferroelectric phase with A = — 3 




2a = b a = b a = 2b 



Figure 13: c-density plots of the antiferroelectric phase with A = — 3 



4 Conclusion 

We demonstrated that local Markov sampling for the 6-vertex model with do- 
main wall boundary conditions indicates that the system develops a macroscopi- 
cal droplet of c- vertices when A < — 1 . For these computations it is not essential 
that the ground state of the 6-vertex model in this phase is doubly degenerate. 
This degeneracy corresponds to the translation by one step in the North-East 
direction on fig. ^] 

This degeneracy is important in the computation of correlation functions 
and other observables. The two ground states correspond to the two parts of 
the graph of states which are connected by a "very narrow neck" in the limit 
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N — ► oo. We will address this problem in the next publication. The existence of 
the droplet can be seen from results of where a different numerical method 
was used. It would be interesting to compare the methods. 

The droplet of c-vertices is similar to the facets in dimer models. The shape 
of the droplet and of the surrounding critical phase is similar to the correspond- 
ing shapes in the dimer model on the square-octagon lattice JT] . 

The local Markov sampling which we used here is equally effective for other 
boundary conditions in the 6-vertex model. Some of the results for more com- 
plicated boundary conditions can be found in 0]. 



A Functions and Implementations 

A.l Main loop 

The following tasks must be completed by the Main loop function: 

1. Import the matrix from a text file 

2. Build flippables list 

3. Set weights 

4. Define p: 



^ max{weight combinations for all flip types} ^ ^ 

5. Loop the following actions, and after a certain defined number of successful 
flips, output a file with the current matrix (and status of the Markov 
Chain) in it. 

(a) Generate a random real, rand, between and 1 

i. If #AiP-"P-™ly + #fli P -down-only + #bi-flip > j tQ 

^vertices — 5 v / 

ii. Otherwise, go to (a). 

(b) Select a random flippable position with probability 

P(selection) — #flip _ up _ only + #flip-down-only + #bi-flip 

by calling the Get Flippable Position function. 

(c) Call Get Weight (which is now scaled by the value of p, to ensure 
that it always returns a value less than 1) to get the probability of 
an up flip and/or a downflip at the flippable location chosen. 

(d) Generate a random real, rand, between and 1. 

i. For up or down-only flips, iff W(Sb) > rand, execute the flip by 
calling the Execute Flip function, else restart main loop. 
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ii. For positions that can flip up and down, iff W(Sb) > rand, 
execute the flip corresponding to Sb, else iff W(S b >) > rand, 
execute the flip corresponding the S b > , else restart main loop. 
In practice, this means that once a vertex which can be flipped 
either way is chosen, simply divide up the probabilities of each 
flip occurring as discussed earlier. 

.2 Execute Flip 

1. If type is high 

(a) Change the entry in the list of Flippables for the vertex chosen to 
make a high flip impossible. 

(b) Define the following positions: 

i. One = the original position = Base 

ii. Two = (+1, +0) = Right 

iii. Three = (+1, +1) = Up Right 

iv. Four = (+0, +1) = Up 

v. Left = (-1, +0) 

vi. Down=(+0, -1) 

vii. UpLcft = (-1, +1) 

viii. UpRight = (+1, +1) 

ix. DownLeft = (-1, -1) 

x. DownRight = (+1, -1) 

(c) Replace 4 parts 

i. Set Contents of Position One = FlipToOne (Position One, High) 

ii. Set Contents of Position Two = FlipToTwo (Position Two, High) 

iii. Set Contents of Position Three = FlipToThrce (Position Three, 
High) 

iv. Set Contents of Position Four = FlipToFour (Position Four, 
High) 

(d) If Up, Down, Right, or Left Positions become flippable (call Get Is 
Flippable on each to check), add them to the Flippables List. 

(c) Call Fix Low End. 

2. If type is low 

(a) Change the entry in the list of Flippables for the entry chosen to 
make a low flip impossible. 

(b) Define the following positions: 

i. One = (-1, -1) = Down Left 

ii. Two = (+0, -1) = Down 
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iii. Three = the original position = Base 

iv. Four = (-1, +0) = Left 

v. Right = (+1, +0) 

vi. Up = (+0, +1) 

vii. UpLcft = (-1, +1) 

viii. UpRight = (+1, +1) 

ix. DownLeft = (-1, -1) 

x. DownRight = (+1, -1) 

(c) Replace 4 parts 

i. Set Contents of Position One = FlipToOne (Position One, Low) 

ii. Set Contents of Position Two = FlipToTwo (Position Two, Low) 

iii. Set Contents of Position Three = FlipToThrcc (Position Three, 
Low) 

iv. Set Contents of Position Four = FlipToFour (Position Four, Low) 

(d) If Up, Down, Right, or Left Positions become flippable (call Get Is 
Flippablc on each to check), add them to the Flippables List. 

(e) Call Fix High End. 

A.2.1 Fix High End 

1. Define the following positions: 

(a) HighCrcateDownLcft = (-1, -1) 

(b) HighDeleteDown = (+0, -1) 

(c) HighDeleteDownDownLeft = (-1, -2) 

(d) HighDeleteDownRighLeft = (-2, -1) 

(e) HighDeleteLeft = (-1, +0) 

2. If HighCrcateDownLcft is flippable, add it to the High Flippables List. 

3. Delete 4 potential flippables on the High End. 

4. If any of HighDeleteDown, HighDeleteDownDownLeft, HighDeleteDown- 
LeftLcft, or HighDeleteLeft exists in the High Flippables list, remove them 
from the list. 

A.2.2 Fix Low End 

1. Define the following positions: 

(a) LowCrcatcUpRight = (+1, +1) 

(b) LowDeleteUp = (+0, +1) 

(c) LowDeleteUpUpRight = (+1, +2) 
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(d) LowDeleteUpRighRight = (+2, +1) 
(c) LowDeleteRight = (+1, +0) 

2. If LowCreateUpRight is flippable, add it to the Low Flippables List. 

3. Delete 4 potential flippables on the Low End. 

4. If any of LowDeletcUp, LowDclcteUpUpRight, LowDcletcUpRightRight, 
or LowDeleteRight exists in the Low Flippables list, remove them from 
the list. 

A.3 Get Weight 

1. Get contents of four surrounding positions if the flip occurred. 

2. Multiply weights together corresponding to the contents of the four posi- 
tions. 

3. Multiply (new weight configuration product) * p 
A.4 Flip To 

FlipTo functions take a position and a type and return what the vertex at 
the given position would be after the flip of the type specified. 

1. FlipToOnc (Position, Type) 

(a) If type is high: 

If vertex was a\, it will be C\; if it was c 2 , it will be a 2 

(b) If type is low: 

If vertex was Ci, it will be a\\ if it was a 2 , it will be c 2 

2. FlipToTwo 

(a) If type is high: 

If vertex was 6 2 , it will be c 2 ; if it was c\, it will be b\ 

(b) If type is low: 

If vertex was c 2 , it will be 6 2 ; if it was b\, it will be c\ 

3. FlipToThree 

(a) If type is high: 

If vertex was a 2 , it will be c\\ if it was c 2 , it will be a\ 

(b) If type is low: 

If vertex was Ci, it will be a 2 ; if it was a\, it will be c 2 

4. FlipToFour 

(a) If type is high: 

If vertex was &i, it will be c 2 ; if it was c\, it will be 6 2 
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(b) If type is low: 

If vertex was c 2 , it will be &i; if it was 6 2 , it will be c\ 

A. 5 Get Flip Position 

1. Generate a random integer between 1 and the total number of flippable 
positions; that is, the number of up-flip only plus the number of down-flip 
only plus the number of bi-flips. 

2. Choose the corresponding element in the Flippable Positions list to the 
random number chosen. 

A.6 Get Is Flippable 

Get Is Flippable should check the status of a position to determine if it is 
flippable. 

1 . High fiippables must be a\ or c 2 vertices and must have empty upper right 
corners (upper right corner must be a 2 or c 2 ). High flippable positions 
must have an x axis coordinate that is less than or equal to the width of 
the matrix 1 (where 0,0 is the origin) and a y axis coordinate that is less 
than or equal to the height of the matrix 1. 

2. Low fiippables must be a\ or c\ vertices and must have empty lower left 
corners (lower left corner must be a 2 or ci). Low flippable positions must 
have an x axis coordinate that is no less than 1 (where 0,0 is the origin), 
and a y axis coordinate that is no less than 1. 
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B Images of the N = 1000 matrix 




Figure 15: N = 1000 c-vertex density plot for the antifcrroclectric phase with 
A = 3 and 2a = b 
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